Create figures for the clinical paper

All data is coming from tables for the paper

# labels for time points
timeLab  = c("Baseline", "Post Run-in", "Week 8")

# set colors
recistCol = c("Non-responder" = "red","Responder" = "blue")
cbrCol = c("CR,PR,SD" = "blue","PD" = "red")
subtypeCol = c("HR+" = "red", "TNBC" = "blue")
timeCol = setNames(c("black", 'grey45', 'grey75'), timeLab)

# set root folder
rootDir = "./"

# Load data patient data
dta_pts = read.csv(paste0(rootDir,"input_tables/table_patient_info.csv"), row.names = 1)
dta_pts$ID = rownames(dta_pts)

#==========================
# file with TILs and PDL1 data
cytof_dat = read.csv(paste0(rootDir,"input_tables/table_cytof_marker_expression.csv"))

#markers = c("CD40", "CD103", "CCR5", "CD141")
markers = colnames(cytof_dat)[4:ncol(cytof_dat)]
#===========================
# add patient info 
cytof_dat = cbind(cytof_dat,  dta_pts[cytof_dat$patientID,])

# take only baseline and time 1 values
df_BvT1 <- cytof_dat %>% filter(timePoint %in% c('Baseline','Post Run-in'))

# take only time 1 and time 2 values
df_T1vT2 <- cytof_dat %>% filter(timePoint %in% c('Week 8','Post Run-in'))

P-values

Wilcoxon test for each time point comparing responders vs non-responders. And paired Wilcoxon p-values for baseline vs post run-in and post run-in vs week 8. Then take differences between all time points and compare responders vs non-responders to see if there is significant difference in changes after treatment

# Calculate p-values from Wilcoxon test: paired for time points and not paired for response
res = c()
# split by clusters
cytof_dat_byClust = split(cytof_dat, f = factor(cytof_dat$cluster))
for(j in names(cytof_dat_byClust))
{
  #
  d = cytof_dat_byClust[[j]]
  for (k in markers)
  {
    # create per patient matrix with three time points in columns
    # some values would be NA
    mat <- data.frame(matrix(nrow = length(unique(d$patientID)),
      ncol = 3, dimnames = list(unique(d$patientID),
                      levels(factor(d$timePoint)))), check.names = F)
    # fill in matrix
    for(i in levels(factor(d$timePoint)))
    {
      d1 = d %>% filter(timePoint == i)
      mat[d1$patientID, i] = d1[,k]
    }
    # add response
    mat$RECIST = cytof_dat[rownames(mat),"RECIST"]
    # create treatment effect 
    # substract baseline from post run-in and week 8 
    mat$delta_T1_B =  mat[,2]-mat[,1]
    mat$delta_T2_B =  mat[,3]-mat[,1]
    # and week 8 from post run-in 
    mat$delta_T2_T1 =  mat[,3]-mat[,2]
    
    # running paired test
    # Baseline vs post run-in
    p1 = wilcox.test(mat[,1], mat[,2], paired = T)$p.value
    # post run-in vs week 8
    p2 = wilcox.test(mat[,2], mat[,3], paired = T)$p.value
    
    # test by response
    # for baseline 
    p3 = wilcox.test(Baseline ~ RECIST, data = mat)$p.value
    # for post run-in 
    p4 = wilcox.test(mat[,2] ~ mat[,"RECIST"], data = mat)$p.value
    # for week 8
    p5 = wilcox.test(mat[,3] ~ mat[,"RECIST"], data = mat)$p.value
    # compare differences in treatment effect in responders vs non-responders
    p6 = wilcox.test(delta_T1_B ~ RECIST, data = mat)$p.value
    p7 = wilcox.test(delta_T2_B ~ RECIST, data = mat)$p.value
    p8 = wilcox.test(delta_T2_T1 ~ RECIST, data = mat)$p.value
   
    # combine all p-values and add means
    res = rbind(res, c(marker = k, cluster = j, BvsT1 = round(p1,3),
                       T1vsT2 = round(p2,3), B_RvsNR = round(p3,3), 
                       T1_RvsNR = round(p4,3), T2_RvsNR = round(p5,3), 
                       d_T1_B_RvsNR = round(p6,3), 
                       d_T2_B_RvsNR = round(p7,3), 
                       d_T2_T1_RvsNR = round(p8,3),
                       mean_d_T1_B_R = round(mean(mat %>% filter(RECIST == "Responder") %>% select(delta_T1_B) %>% unlist(), na.rm = T),3),
                       mean_d_T1_B_NR = round(mean(mat %>% filter(RECIST == "Non-responder") %>% select(delta_T1_B) %>% unlist(), na.rm = T),3),
                       mean_d_T2_B_R = round(mean(mat %>% filter(RECIST == "Responder") %>% select(delta_T2_B) %>% unlist(), na.rm = T),3),
                       mean_d_T2_B_NR = round(mean(mat %>% filter(RECIST == "Non-responder") %>% select(delta_T2_B) %>% unlist(), na.rm = T),3),
                       mean_d_T2_T1_R = round(mean(mat %>% filter(RECIST == "Responder") %>% select(delta_T2_T1) %>% unlist(), na.rm = T),3),
                       mean_d_T2_T1_NR = round(mean(mat %>% filter(RECIST == "Non-responder") %>% select(delta_T2_T1) %>% unlist(), na.rm = T),3)))
  }
}
res = res[order(res[,'marker']),]

datatable(res)
# make boxplots for markers/clusters that have p-value < 0.05 in d_T1_B_RvsNR
sets = data.frame(res)
sets = sets %>% filter(d_T1_B_RvsNR < 0.05 & cluster%in% c('B','CD4_T','CD8_T','DNT','DPT','NK_like'))

for (i in 1:nrow(sets))
{
  k = sets[i,'cluster']
  j = sets[i,'marker']
d = cytof_dat %>% filter(cluster == k)
d = CreateAllFacet(d,"RECIST")
  df = cbind(protein = d[,j], 
             d[,c("patientID", "timePoint", "RECIST","Subtype", "facet") ] )

    p = ggpaired(df , x="timePoint",y="protein", id = "patientID", 
                 line.size = 0.4, title = paste0(j,". ", k), 
                 line.color = "RECIST",palette =  recistCol, xlab="Time") +
      facet_wrap("facet") +
      scale_x_discrete(labels=timeLab)
    
    print(p)
}

Baseline vs Post Run-in

#ggpaired(df_BvT1 %>% filter(cluster %in% c("GMDSC")),
for (i in markers)
{
  p = ggpaired(df_BvT1, x="timePoint",y=i, id = "patientID", line.color = "RECIST", line.size = 1,
   palette =  recistCol, title = paste0(i,". Baseline vs Post Run-in"), xlab="Time", ylab="MMI") + 
    scale_x_discrete(labels=timeLab)+ 
    facet_wrap("cluster", scales = 'free') +
    theme(axis.text.x=element_text(angle = 45, hjust = 1))

  print(p)
}

Post Run-in vs Week 8

for (i in markers)
{
  p = ggpaired(df_T1vT2, x="timePoint",y=i, id = "patientID", line.color = "RECIST", line.size = 1,
   palette =  recistCol, title = paste0(i,". Post Run-in vs Week 8"), xlab="Time", ylab="MMI") +
#    scale_x_discrete(labels=timeLab)+
    facet_wrap("cluster", scales = 'free') +
    theme(axis.text.x=element_text(angle = 45, hjust = 1))

  print(p)
}

sessionInfo()
## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] scales_1.2.1      gridExtra_2.3     DT_0.26           data.table_1.14.6
## [5] ggrepel_0.9.2     ggpubr_0.5.0      ggplot2_3.4.0     dplyr_1.0.10     
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.0  xfun_0.36         bslib_0.4.2       purrr_1.0.1      
##  [5] carData_3.0-5     colorspace_2.0-3  vctrs_0.5.1       generics_0.1.3   
##  [9] htmltools_0.5.4   yaml_2.3.6        utf8_1.2.2        rlang_1.0.6      
## [13] jquerylib_0.1.4   pillar_1.8.1      glue_1.6.2        withr_2.5.0      
## [17] DBI_1.1.3         lifecycle_1.0.3   stringr_1.5.0     munsell_0.5.0    
## [21] ggsignif_0.6.4    gtable_0.3.1      htmlwidgets_1.6.1 evaluate_0.19    
## [25] labeling_0.4.2    knitr_1.41        fastmap_1.1.0     crosstalk_1.2.0  
## [29] fansi_1.0.3       highr_0.10        broom_1.0.2       Rcpp_1.0.9       
## [33] backports_1.4.1   cachem_1.0.6      jsonlite_1.8.4    abind_1.4-5      
## [37] farver_2.1.1      digest_0.6.31     stringi_1.7.12    rstatix_0.7.1    
## [41] grid_4.2.2        cli_3.4.1         tools_4.2.2       magrittr_2.0.3   
## [45] sass_0.4.4        tibble_3.1.8      tidyr_1.2.1       car_3.1-1        
## [49] pkgconfig_2.0.3   ellipsis_0.3.2    assertthat_0.2.1  rmarkdown_2.19   
## [53] rstudioapi_0.14   R6_2.5.1          compiler_4.2.2